On the high dimensional Bak-Sneppen model 
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We report on extensive numerical simulations on the Bak-Sneppen model in high dimensions. We 
uncover a very rich behavior as a function of dimensionality. For d > 2 the avalanche cluster becomes 
fractal and for d > 4 the process becomes transient. Finally the exponents reach their mean field 
values for d = d c = 8, which is then the upper critical dimension of the Bak Sneppen model. 
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The Bak-Sneppen (BS) model 0, since its introduc- 
tion, has attracted much attention in statistical physics. 
Thanks to its simplicity, it has lead to a much deeper 
understanding of the nature of self organized criticality 
|^| and of extremal dynamics H in general. The rules of 
the BS dynamics, in d dimensions, are very simple: the 
state of the model is completely denned by L d numbers 
fi arranged on a d-dimensional lattice of edge size L. At 
every time step the smallest of these numbers and its 
2d nearest neighbors are replaced with new uncorrelated 
random numbers, drawn from the uniform distribution. 
Such very simple dynamics, based on the selection of the 
global minimum, is generally called extremal dynamics. 
It was first introduced in invasion percolation Q and it 
results in a remarkably rich and interesting critical be- 
havior. 

The self-organized critical nature of the BS model (as 
well as of other extremal models) is revealed in its abil- 
ity to naturally evolve towards a critical state where al- 
most all of the variables fi are above a threshold f c . 
The dynamics in this state is characterized by scale- 
free bursts of activity or avalanches, which form a hi- 
erarchical structure of sub- avalanches within big- 
ger avalanches. This critical state is described by crit- 
ical exponents. The distribution of avalanche duration 
s behaves as a power law P(s) ~ s~ T with exponent 
r. An avalanche of duration s covers a number of sites 
V(s) ~ s M . The set of avalanche sites is generally frac- 
tal V(s) ~ R(s) Df , where R(s) is the gyration radius 
and Df is the (spatial) fractal dimension. The active 
site (the one with the global minimum fi) has a dynam- 
ical wandering that can be described in terms of return 
times: The distribution Pf(t) ~ t~ Tf of first return times 
is characterized by an exponent Tf , whereas the distribu- 
tion of all return times P a {t) ~ t~ Ta defines the exponent 
r a . For a random walk these two exponents take the val- 
ues Tf = 3/2 and r a = 1/2 in d = 1, whereas for d > 2 
one finds 77 = r a = d/2. Note that with respect to pre- 
vious literature 0,0], our notation is slightly different: 
The exponent fi in ref. Q is defined as fi = dj D, where 
D was called "fractal dimension". Here we regard /j as 
a fundamental exponent (not a composite one) and re- 



serve the name fractal dimension to describe the spatial 
geometrical properties of the avalanches. This choice, as 
we shall see, avoids ambiguities which can arise in higher 
dimensions, where the avalanche cluster becomes fractal. 

A scaling theory was proposed in which shows that 
scaling relations allow to reduce the number of critical 
exponents to two, for example fi and r. Numerical sim- 
ulations in d = 1 and 2 fully confirms the validity of the 
scaling theory The mean field limit, formally cor- 
responding to the limit d — > 00, has also been solved 
exactly [B: the exponents, in this limit take the values 
t = Tf = r a = 3/2, [i — 1 and Df — 4. Finally, it was 
recently shown [^|J^] that a further non-trivial relation, of 
a different nature, exists between \x and t. This suggests 
that the BS universality class is characterized by a single 
exponent, e.g. ^(d), as a function of d. 

In this Letter we analyze the behavior of the BS model 
as a function of dimensionality. The understanding of 
the critical behavior of a model, as a function of di- 
mensionality, is a central issue in statistical physics. In 
particular the identification of the upper critical dimen- 
sion d c , above which the mean field picture applies, is 
of great importance. Indeed it allows to understand the 
behavior of a finite dimensional system using the pow- 
erful tools of dimensional (e) expansion. In equilibrium 
statistical mechanics, this is almost routine work, but for 
non-equilibrium systems it is a challenging issue. For this 
reason the understanding of the behavior of simple mod- 
els as a function of dimensionality is of great importance. 

We present extensive numerical simulations for the BS 
model which show that i) the upper critical dimension is 
d c = 8 ii) for d > 2 the avalanches are no more compact 
Df < d Hi) in d — 3 we find fi < r a < 1 and iv) for 
d > 4 the process becomes transient i.e. r n = r/ > 1. 
The BS model then shows a quite rich behavior, with 
four qualitatively different regimes (d < 2, 2 < d < 4, 
4 < d < 8 and d > 8), as a function of dimensionality. 
In particular we find that the relation /1 = r a holds only 
up to d — 2 (note that this relation is in any case prob- 
lematic close to the mean field limit \i = 1, r a = 3/2). A 
possible interpretation of our results is that "geometric" 
exponents, such as r a , Tf and Df become independent 
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from "avalanche" or "memory" (see later) exponents for 
d > 2 and the number of independent exponents changes 
with dimension. We shall first discuss in detail the nu- 
merical procedure, defining operationally the quantities 
we measure. Then we present the numerical results and 
finally we discuss them. 

The BS dynamics can be simulated very efficiently us- 
ing a tree like search and replace algorithm B. This 
decreases greatly computation times, so that memory is 
the only limitation to the system sizes studied. The nu- 
merical procedure was first tested in d — 1 and d = 2 and 
we found complete agreement with refs. ||||. In order 
to compute the exponents fi and r, following ref. [flo|| , 
we introduce an age variable ki on each site. The age 
ki of site i measures the time elapsed since the last up- 
date of the variable This method enables us to give 
a precise evaluation of the exponent /x. Indeed at each 
time the sites with age ki less than s identify the current 
avalanche of duration s. Therefore V(s) is simply ob- 
tained counting the sites with ki < s. Age variables also 
allow for a determination of the avalanche exponent r. In 
systems evolving by an extremal dynamics, it has been 
found |5 10 that the probability that the global minimum 
fi occurs on a site with age ki behaves as p(ki) ~ k^ a . 
The exponent a > implies that older sites are less likely 
to be selected. This reflects the fact that a very old site 
has survived many selections, and therefore it has very 
likely a high value of the variable fi. The more it sur- 
vives the higher is its fi, possibly even greater than the 
threshold / c , in which case the site will never be selected 
before a new update (occurring if one of its neighbors is 
selected) . The exponent a is related to the avalanche ex- 
ponent by a — 3 — r (|llj. The advantage of measuring a, 
with respect to a direct measure of r, is that the statis- 
tics of the former is much richer than that of the latter in 
the same simulation. We measured r in both ways and 
found good agreement (which also supports the validity 
of the relation a = 3 — r) . Since statistical uncertainty of 
the exponent a is much less than that of r we report here 
only the value r = 3 — a. The return times exponents are 
measured in the usual way @: let tf \ k = 0,1,... be 
the (return) times when site i is visited (t^ < t[ k ). 
The first return exponent is obtained from the statistics 
of _ ft 1 whereas the all returns exponent is ob- 



tained from tj*^ — t\"'. Finally, in order to obtain the 
fractal dimension Df, we compute the gyration radius of 
avalanches of size V(s). 

The numerical results are summarized in Table |. As 
shown in Fig.[l], our numerical results for /i and r are in 
good agreement with the exact relation recently found in 
ref. 0. This is an important consistency check for the 
reliability of the simulations. In d = 2 we find a slight 
difference between r a and /x which is however within er- 
rorbars. The data are also plotted in Fig.|| for complete- 
ness. For d = 3 the difference between r Q and fj, is much 
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larger than the error bars (see Fig.||). In Fig.|| we plot 
the exponents as a function of dimensionality. For d < 3 
the process is recurrent, since r Q < 1, i.e. each visited 
site, in an infinite system, is visited again an infinite 
number of times (or, stated differently, a selected site 
will be selected again with probability one). Instead for 
d > 4 the process becomes transient i.e. r a = t/ > 1: 
each site, in an infinite system, is visited a finite number 
of times (there is a finite, smaller than one, probability 
that a selected site will be selected again). The relation 
r a + Tf = 2 if r a < 1 and r a = r/ if r a > 1 (a classical 
result from renewal theory |l2|]) is always respected, the 
former holding for d < 3. The return time statistics, in 
the BS model, is determined both by memory effects and 
by geometry. Activity can return to site i either because 
all sites with variables larger than fi are eliminated by 
the BS dynamics (which we call a memory return) or 
because the activity returns close to site i and fi — > f[ 
is updated (geometric return). As the dimension d in- 
creases, geometric returns become less and less relevant 
(see later) and for very large dimensions one expects to 
recover the mean field result r a = rt =3/2. We find that 
the mean field limit holds for d > d c = 8. 

We now turn to the discussion of the onset of fractal- 
ity of avalanche clusters for d > 2. At d = 2, we found 
numerically that the fractal dimension is slightly smaller 
than Df = 2. We argue, however, that small size effects 
occur and that Df = 2 holds. Indeed, as the system size 
L is increased the slope of the curve \ogV(s) vs log i?(s) 
increases, suggesting that avalanches are compact. The 
occurrence of small size effects can be understood analyz- 
ing the dynamics of the growth of the avalanche cluster 
in d dimensions. The growth probability distribution of 
the cluster is characterized by memory. The problem of 
cluster growth with memory was recently investigated in 
ref. 11311 . In a standard model of growth, new sites can 
be added only at the surface of the cluster. In this case, 
it was shown that when a > 1 the cluster has a fractal 
dimension Df < d [jlji. The difference between cluster 
growth with memory and the avalanche cluster growth 
in the BS is that in the latter the activity can return to 
bulk sites. At the initial stages of growth of BS clusters, 
however, the analogy is useful, since the ratio between 
the number of surface sites and that of bulk sites is ~ 2d 
so that returns to the bulk are rare. In d = 2, as more 
and more returns to the bulk occur the surface to volume 
ratio decreases and the cluster becomes compact (for a 
fractal cluster the ratio of perimeter sites to bulk sites is 
finite). Therefore, while small s avalanches look fractal, 
for larger sizes they become more and more compact. 
These small size effects are also visible in the distribu- 
tion of first return times (see also Q). In d > 3 the 
same size effect should appear, however in this case we 
find a fractal dimension which is definitely smaller than 
d p3|. Returns to the bulk are fewer in d > 3 than in 
d = 2 and more importantly the topology of a fractal 
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cluster is very different: while the 2-dimensional cluster 
is characterized by a distribution of holes of all sizes, a 
fractal in d > 3 is most likely a ramified object. Returns 
to the bulk can fill the holes of the d = 2 cluster but it is 
much more difficult for them to turn the branched frac- 
tal structure of the d > 3 avalanche into a compact one. 
For this reason, we believe that our data are compatible 
with the occurrence of a compact cluster in d = 2 and 
with Df < d when d > 3. In any case, numerical data 
are relatively stable with changes of the system size and 
we could not detect any small size effect such as the one 
discussed above for d = 2. 

Geometric returns arise because the avalanche cluster 
has many self-intersections: indeed the self-intersection 
set has a fractal dimension Dj = 2Df — d > for d < 8 
|fl5f . Note that Dj is smaller than D f which means that 
the larger the avalanche the smaller is the fraction of in- 
tersection sites. Moreover, as mentioned earlier, geomet- 
ric returns become less and less relevant as d increases 
since Df — Di = d — Df also increases with d. Above 
d = d c = 8 the fractal dimension Df = 4 attains its 
mean field value and the avalanche has no self intersec- 
tions (Di < 0). 

The fractality of avalanches implies that some scaling 
relations must be modified. For example the jump prob- 
ability distribution p(r) was analyzed in ref. 0] under the 
hypothesis of a compact cluster. It is defined as the prob- 
ability that the activity jumps in one time step to a site 
at a distance r from the current site and it falls off with 
a power law behavior: p(r) ~ r _7r . The exponent tt was 
related to r and D in ref. Q by tt = 1 + D(2 — r), under 
the hypothesis of a compact cluster. The same argument 
as in ref. led us to the expression 

Df(2-T) 

ir = l + ^ '-. (1) 

M 

It is interesting to observe that, for the values quoted 
in Table |[ tt is always slightly larger than 3, so that the 
second moment of p(r) is finite, and its mean field value is 
tt = 3. Were it not for the correlations, the activity would 
perform a random walk and it would be transient already 
for d > 2. Correlations induced by memory effects, shifts 
the onset of transient behavior to d > 4. 

In order to check the full consistency of the results 
shown in Table |, we are also exploring an alternative 
way to move away from the mean field limit [^6[ . Within 
a d = 1 system, we choose the "nearest neighbors" of 
the active site at random over the lattice, with a prob- 
ability which is a power law decreasing function of the 
distance from the active site, with exponent to. Prelim- 
inary results show that the mean field limit is recovered 
when u> — * 1. Varying lj > 1 again three different crit- 
ical regimes appear. In particular we find a region of 
the values of u> where activity is recurrent (r a < 1) but 
p =/= r a and a region where activity is transient (r Q > 1) 
but p < 1. Moreover we find that, whenever p takes the 



same values listed in Table Q, also the other quantities 
take on the same corresponding values (apart of course 
from Df which is always smaller than 1). Finally, work- 
ing in d — 1 allows us to simulate extremely large sizes 
and to rule out any finite size effects, thus strenghtcning 
the reliability of the results shown in Table ffl. 

It has been argued that the mean field limit of BS 
can be described as a branching process ||. The value 
d c = 8, which coincides with the upper critical dimension 
of branched polymers, is consistent with this picture |l7[ . 
We believe that close to d = 8 branched polymers give a 
reasonable description of the geometry of the BS process. 
An e = d c — d expansion for the BS model could therefore 
be feasible, also using the recent expansion around the 
mean field solution of ref. J?]]. In this respect, our results 
allow for a prediction of the first coefficient p ~ 1 — 0.017e 
of the e expansion. 

In summary, we have presented numerical results for 
the BS model in high dimensions. These allow to con- 
clude that for d > 2 avalanches become fractal and for 
d > 4 the process becomes transient. Finally the mean 
field limit is reached at the upper critical dimension 
d c = 8. This behavior, with three different non-trivial 
regimes of criticality, is substantially richer that that of 
equilibrium statistical models, where non-trivial behav- 
ior occurs only in one region of dimensionality, limited by 
the lower and the upper critical dimension. Our results 
also call for a close analysis of the scaling theory of ref. 

, which has been developed under the implicit assump- 
tions of a recurrent process with compact avalanches. 
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TABLE I. Exponents of the BS model for different dimen- 
sions d. An upper bound on the error is ±0.01 for safety, even 
if for some of the listed exponents confidence is greater. In 
particular, for d — 7, the values are distingushable from the 
d — 8 ones. The last line gives the size L — 2 n of the edge of 
the hypercube used for the simulations. 



FIG. 1. Avalanche exponent r vs. /i for different dimen- 
sions (crosses). The continuous line is the exact relation be- 
tween the two exponents as from ref.[7]. The data from the 
simulations and the exact relation are in excellent agreement. 

FIG. 2. Log-log plot of the number of covered sites V(t) 
(dashed line) and the (inverse) all return times distribution 
P a {t) (full line) as a function of time t. The results in a) d — 3 
show that the slope of the two quantites are clearly different 
(hence n =fc r a ) b) d = 2 show instead full compatibility of the 
slopes with fi = r a . 

FIG. 3. The exponents fi, r, r a and r/ as a function of 
dimensionality. Dashed lines at 1 and 1.5 have been drawn 
for reference. 
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